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SUMMARY 

A nonlinear finite element model of a rotor bearing system was linearized at some 
instant of time. The stiffness and damping coefficients were assumed to be constant 
over a short time period at the instantaneous values. A numerical stability analysis 
was applied to the linearized model to determine the stability limits of the forward, 
backward, and centered Euler; the fourth-order Runge-Kutta; the Milne; and the Adams 
methods of numerical integration of the set of equations of motion. 

The stability analysis obtained the error growth rate (the ratio of the amplitudes of 
the errors for two succeeding time steps) as a function of a nondimensional time step 
(the product of the time step and the natural frequency) and of the damping ratio for each 
mode. The absolute error growth rate was independent of the actual transient; the rela- 
tive error growth rate was not. The relative error growth rate was normalized by the 
transient solution of the perfectly balanced rotor. 

For any mode with a nonnegative damping ratio, the absolute error growth rate 
must be less than or equal to 1. (If not, the numbers used in the calculation will be- 
come too large for the computer.) If a mode is to be calculated accurately, the relative 
error growth rate must also be less than or equal to 1. The lower frequency modes 
are the only modes that are physically meaningful, although all the modes are inher- 
ently present in the calculation of any transient. 

Usually, the hipest frequency mode possible in the finite element model determines 
the threshold of numerical stability. Therefore, the number of mass elements in the 
finite element model should be minimized. Increasing the damping of a mode can cause 
it to become numerically unstable. An example considered is that of a uniform shaft on 
rigid bearings, with 10 mass elements, and operating at approximately the first critical 
speed. The maximum time step for the Runge-Kutta, Milne, and Adams methods is 
that which corresponds to approximately 1 degree of shaft movement. To calculate the 
shaft motion to several revolutions would thus require thousands of time steps. 



INTRODUCTION 


The simulation of transient rotor dynamics by a computer is done by one of two 
methods. The first is the modal method. This method is applicable to linear prob- 
lems. It consists of integrating the equations of motion for each mode separately and 
then summing them to get the final result. The second method is the finite element 
method. This method is applicable to nonlinear problems. It consists of modeling the 
rotor dynamics system by a finite number of elements and then integrating the equations 
of motion for all of these elements simultaneously. 

References 1 and 2 use the finite element approach. Reference 1 is the more gen- 
eral of the two codes. Besides simulating more features of the rotor dynamics sys- 
tem, it includes several choices of numerical method for integrating the equations of 
motion. Reference 1 concluded that the Adams- Moulton predictor- corrector integration 
technique gave the best combination of computational speed and accuracy. Application 
of this code to other problems, such as that of reference 3, has been frustrated by nu- 
merical stability problems. 

Reference 2 is the more restrictive of the two codes. It uses a modified (centered) 
Euler technique to integrate the equations of motion. This code has been applied suc- 
cessfully to the problem of reference 3. The question arises as to why the code of ref- 
erence 2 worked and that of reference 1 did not. A partial answer to this question can 
be found in the numerical stability of the various integration techniques. 

This report discusses the numerical stability of some typical integration techniques 
as applied to transient rotor dynamics. The set of equations simulating the transient 
rotor dynamics would be classified mathematically as a stiff Set of equations. (A set of 
equations is stiff if the ratio of the largest to the smallest eigenvalue is large, i.e. , 
more than 100 to 1.) Reference 4 discusses numerical stability and proposes a numeri- 
cal integration technique, Gear's method, for a general set of stiff equations. This 
method was applied to the code of reference 1 but required too much computing time. 

The numerical stability analyses of references 4 and 5 were applied to some typical 
numerical methods presented in reference 6. The present study examined the numeri- 
cal stability of the forward, backward, and centered Euler methods; the fourth-order 
Runge-Kutta method; the Milne method; and the Adams method as applied to a transient 
rotor dynamic simulation. 

The finite element method divides the shaft into a number of axial elements. The 
acceleration of the elements is related to the sum of the forces on the elements. These 
forces are basically the elastic force, the drag force, and the inertial unbalance force. 
The elastic force is related to the shaft stiffness and therefore to the displacements of 
the elements. The drag force is related to the velocity of the elements. The unbalance 
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force is independent of the displacement and velocity of the elements and is therefore 
the "forcing function." 

The numerical integration technique replaces the time- dependent set of differential 
equations of motion by a set of difference equations. When the time step is small, the 
solution to the difference equations approximates the solution to the differential equa- 
tions. When the time step in the difference equations exceeds a critical value, the so- 
lution no longer approximates the solution of the differential equations of motion. The 
relative error between the two solutions increases with time. The numerical integra- 
tion technique is unstable for time steps greater than the critical value since a small 
error will increase to a large error. 

This critical time step is not generally known since the solution to the differential 
equations of motion is not known. If the problem is linearized in time (i. e. , the stiff- 
ness and damping of the elements are assumed to be constant at the instantaneous val- 
ues), a modal analysis can be applied to the instantaneous mode shapes and frequencies 
for both the differential and difference equations. The relative error in the amplitude 
of each mode must not increase with time. The stability analysis applies for the lin- 
earized problem; therefore, an extension must be made between the linearized and non- 
linear problems (i, e. , if the linearized problem is unstable, the nonlinear problem is 
also unstable). 

Since the problem is linearized, the absolute error can be calculated independently 
from the actual simulation. The error must satisfy the same set of homogeneous equa- 
tions, but it has a different forcing function (the rate at which the computer generates 
errors). If the computer rounds the last significant figure (rather than truncating), the 
forcing function, on the average, is zero. The error analysis then reduces to an initial- 
value problem that can be solved analytically for each numerical integration method. 

In order to determine the relative error, the absolute error must be normalized by 
the solution to the linearized nonhomogeneous differential equations of motion. This 
solution depends on the specific transient. In order to generalize the results, the most 
conservative transient is used to normalize the absolute error; that is the solution to 
the homogeneous differential equations of motion. The solution to the homogeneous 
equation of motion (i.e. , the shaft in perfect balance) will yield the smallest vibrations 
and therefore the largest relative error. 


ANALYSIS 

This analysis assumes a model of a rotor bearing system that is linearized at some 
instant of time and neglects both torsional windup and gyroscopic moments. Figure 1 
shows a model of the shaft with n finite elements. The complex number representation 
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of the radial displacements of the shaft centerline from the axis of rotation for the i*^^ 
finite element is r. (Symbols are defined in the appendix. ) 

If R is a column matrix of the displacements of ttie shaft centerline, if the column 
matrix V is the time derivative of R, and if the column matrix A is the time deriva- 
tive of V, the second^order matrix equation of motion for the rotor bearing system can 
be written as 


MA + CV + KR = F 


where M is a diagonal matrix of the masses, C is a square matrix of the damping 
coefficients, K is a square matrix of the stiffness coefficients, and F is a column ma- 
trix of the forces. The following partitioned matrices can be defined: 

Z = 




Q 



0 

F 


Then the first-order matrix equation of motion becomes U = PZ + Q, where U is the 
time derivative of Z. 

If Z. is defined to be the solution of the eigenvalue equation PZj = , the solu- 

tion of the homogeneous differential equation of motion can be expressed as a modal 
series: 


Z = 




z. 

j 
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where aj, the complex amplitude of the j mode, is 


a^ = aj(0)exp(Ajt) 


o 


From the definition of P, 


and 




A.R. 
J J 


-M“^(KRj. + CV^) = XjV^. 


Substituting for Vj yields the damped eigenvalue equation 



M +X.C +K 
J 


R. = 0 
J 


with 2n eigenvalues. If this equation is premultiplied by Rj (Rj conjugate trans- 
posed), it becomes 

X^m. + X.c. + k. = 0 
J J J J J 

where the scalar modal mass, modal damping, and modal stiffness are defined as 


m. 

J 


R*MR. 

= _J 1 

* 

R. R. 

J I 


c. 

J 


* 

R. CR. 

= ^! 1 

* 

R. R. 

J J 


k. 

J 


R.KR. 

-J 1 

* 

R. R. 

J J 


If the natural frequency is defined as 




and the damping ratio as 


c. 



the eigenvalue is 



The frequency is defined to be greater than zero. If the damping ratio is greater 
than zero, the mode is damped; if the damping ratio is less than zero, the mode is am- 
plified. 

The firsl^order matrix equation of motion defines U as a function of Z and t. 

The various numerical solutions of this set of differential equations are solutions of the 
following difference equations: 

Forward Euler method: 

Z(t +h) = Z(t) +hU[Z(t), t] 

Backward Euler method; 

Z(t + h) = Z(t) + hU[Z(t + h), t + h] 

Centered Euler method: 

Z(t +h) = Z(t) +-^|u[Z(t),t] + U[Z(t +h),t +h)]| 
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Runge-Kutta method (fourth order): 


Z(t +h) = Z(t) 


. +2gg+g^ 

6 


where 


= hU[Z(t), t] 


^2 


^3 


= hU 


= hU 


% h 

Z(t) +— , t +- 

2 2 


®2 h 

Z(t) +— , t +- 

2 2 


g4 = hU[Z(t) +g3, t +h] 


Milne method: 


Z(t +h) = Z(t- h) 


+-|u[Z(t +h) 


t + h] 


+ 4U[Z(t),t] +U[Z(t-h),t- 



Adams method: 


Z(t +h) = Z(t) +-lL|9U[Z(t +h),t +h] + 19U[Z(t),t] - 5U[Z(t- h),t- h] 


+ u[Z(t- 2h),t- 2h]j 


From the definition of U, these difference equations can be generalized as 


Z(t+h) = ^ D^Z(t- Zh) +B^Q(t- Ih) 

I 

where I indicates the various terms used in the difference equations and and 
are matrices that are different functions of P and h for each difference equation. The 
exact solution to these difference equations is Z. The computation actually calculates 
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Z(t +h) +E(t +li) |D^[Z(t- Ih) + E(t- Ih)] +B^Q(t- Zh)J + G(t 


h) 


where E is a column matrix of the roundoff error and G is a column matrix of the 
rate of generation of the roundoff error. If the specific computer rounds the last sig- 
nificant figure (rather than truncating), the rate of generation of roundoff error, on the 
average, is zero. Subtracting the exact difference equation from the computational dif- 
ference equation yields the following homogeneous difference equation for the roundoff 
error: 


E(t +h) = ^ D^E(t- ih) 


The matrix can be expressed as a polynomial function of hP 




dk;(hP)‘ 


where dj^^ is a different scalar constant for each integration technique and k indicates 
the various powers of the polynomial. Table I is a tabulation of for the various 
techniques. If E is expressed as a sum of eigenvectors (modal expansion), E becomes 


E = 




fVi 

where e^. is the complex amplitude of the roundoff error of the j mode. Substituting 
this into the difference equation for the roundoff error and noting that Z^ is a linear in- 
dependent eigenvector of P result in the following equation for the amplitude of the 
roundoff error of the j mode: 


e. 

J 




I k 


dk^ (hAj)^ej (t - ^h) 


This homogeneous difference equation has a solution of the form 


ej(t) = e^(0)exp(^jt) 
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where juj is a constant, 
fined as 


If Wj, 


the absolute error growth rate of the 


.th 

J 


mode, is de- 


e. (t + h) 
w. =— ■! 


Wj becomes 


Wj = expOijh) 


If 


Sj, a complex nondimensional time step for the j 


• th 


mode, is defined as s. = hX., 

j J 

this form of solution converts the difference equation into an algebraic equation in 


terms of w^ 


and 


s.: 

J 


A .,-1 

I k 


w. = ) > d. , s. W- 

J ki J J 


fVi 

If the relative error of the j mode is defined as 

e.(t) 
e (t) =-i_ 

J aj(t) 

the relative error growth rate of the mode is 

e (t + h) 
u. = — J 

Substituting the values of a. and e. and using the definitions of s. and w. yield 

3 J 3 3 

Uj = Wj exp(-s) 

The algebraic equation in terms of w^ and Sj is transformed into an equation in terms 

of u. and s. 

3 3 


I k 


rl 
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Both equations are polynomials in terms of either Wj or m with complex coeffi- 
cients that are functions of s. These polynomials can be solved either analytically or 


numerically for Wj 

X., s. becomes 
J 3 


or 


u. 

J 


s. 

in terms of 


s. From the definition of 




and the value of 


s.=hcu. 



The largest absolute value of either or m was used with a numerical contour 
plotting routine in terms of hco^ and ty These contours are shown in figures 2 and 3. 


DISCUSSION 

A nonlinear finite element model of a rotor bearing system can be linearized at any 
instant of time. The stiffness and damping coefficients are assumed to be constant over 
a short time period at the instantaneous values. A damped critical- speed analysis can 
then be applied to obtain the instantaneous mode shapes, natural frequencies, and damp- 
ing ratios. The analysis only applies for the time period for which the linearization is 
valid. For other time periods the analysis must be repeated with "new” stiffness and 
damping coefficients. 

A numerical stability analysis was applied to the linearized finite element model of 
the rotor bearing system to determine the stability limits of the forward, backward, and 
centered Euler; the fourth-order Runge-Kutta; the Milne; and the Adams methods of 
numerical integration of a set of differential equations. The stability analysis obtained 
the error growth rate (the ratio of the amplitudes of the errors for two succeeding time 
steps) as a function of a nondimensional time step (the product of the time step and the 
natural frequency) and of the damping ratio for each mode. Figure 2 shows contour 
lines of the absolute error growth rate and figure 3 shows contour lines of the relative 
error growth rate. The relative error growth rate was normalized by the transient so- 
lution of the perfectly balanced rotor. 

In the finite element model of the rotor bearing system, the higher frequency modes 
are not physically meaningful but are inherently present. The accuracy of the calcula- 
tions in these modes is not important. It is only important that, for a nonnegative 
damping ratio, the absolute error be bounded (i. e. , the absolute error growth rate must 
be less than or equal to 1). If the absolute error growth rate is greater than 1, even- 
tually the numbers used in the calculation will become too large for the computer and a 
computer overflow will occur. This is not too serious for the negative damping ratio 
case (amplification) since the rotor vibrations will be unbounded and the actual displace- 
ments will eventually cause the computer to overflow. 
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The lower frequency modes are the only modes that are physically meaningful. The 
accuracy of the calculations in these modes is important. These modes must have a 
relative error that is bounded (i. e. , the relative error growth rate must be less than or 
equal to 1). If the relative error growth rate is greater than 1, eventually the numbers 
used in the calculation will become meaningless. The relative error growth rate shown 
in figure 3 is for a rotor that is in perfect balance. This will yield the smallest vibra- 
tions and therefore the largest relative error. A rotor that is out of balance will have 
larger vibrations and therefore a smaller relative error. K ttie integration technique is 
to be universal (i. e. , applying to all transients), relative error growth rate as shown in 
figure 3 must be less than or equal to 1. 

In general, increasing the damping ratio while keeping the product of the time step 
and the natural frequency fixed causes the relative error to grow. This is also true for 
the absolute error, with the exception of the centered and backward Euler methods. 

For damping ratios between -0. 5 and 0. 7 for the forward Euler method, between 0. 5 
and oo for the backward Euler method, between 0 and 0.8 for the centered Euler method, 
between 1.7 and °o for the Runge-Kutta method, between 0 and °o for the Milne method, 
and between -oo and -1.7 for the Adams method, the relative error growth rate is 
greater than 1. 


EXAMPLE 

As an example, consider a uniform shaft on rigid bearings. The damping ratio is 
zero for all modes; thus, the absolute and relative error growth rates are equal. The 
natural frequency for a continuous model of this system is 

_ .2 

CU. = 1 CJ, 

J ^ 1 

As an order-of-magnitude approximation, the natural frequencies of the finite element 
model are approximately equal to the natural frequencies of the continuous model. If 
there are n elements in the model, the hipest natural frequency is 

The n‘'“ mode (not usually desired from the calculations but inherently present) 
determines the threshold of the numerical instability. From figure 3 the forward Euler 
method is always unstable. The backward and centered Euler methods are never un- 
stable. The Runge-Kutta method is unstable for a time step greater than approximately 
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3 


h > 

2 

n 

The Milne method is unstable for a time step greater than approximately 


2 

n 

The Adams method is unstable for a time step greater than approximately 


2 

n 

If the shaft is running at approximately the first natural frequency and if there are 
10 elements, the calculation must be done at least every 1.5, 1, and 0. 5 degrees of 
shaft rotation for the Runge-Kutta, Milne, and Adams methods, respectively. To cal- 
culate the shaft motion to several revolutions would require thousands of time steps. 

If these time- step limits are exceeded by a factor of 2, the error growth rates are 

approximately 10, 3, and 1. 5 for the Runge-Kutta, Milne, and Adams methods, respec- 

—8 

tively. This means that a relative error of 10 will grow to a value of 1 in less than 
8, 17, and 46 time steps for each of these methods, respectively. 

The preceding example illustrates the fact that great care should be taken to mini- 
mize the number of elements in the finite element modal, thus maximizing the time 
step. This is important not only to save calculation time but also because the time step 
has a lower limit dictated by the roundoff error. 


CONCLUSIONS 

The numerical stability of the Euler, Runge-Kutta, Milne, and Adams integration 
techniques as applied to transient rotor dynamics was analyzed. The following conclu- 
sions were drawn: 

1. The highest frequency mode possible in the finite element modal usually deter- 
mines the threshold of numerical stability. Therefore, the number of mass elements 
in the finite element model should be minimized. 

2. Increasing the damping ratio of a mode can cause it to become numerically un- 
stable. 
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3. Numerical stability for the Runge-Kutta, Milne, and Adams methods for a uni- 
form shaft on rigid bearings with 10 mass elements (operating approximately at the first 
critical speech is to make a calculation at less than 1 degree of shaft rotation. 

In general, the product of the time step and the natural frequency, for all modes pos- 
sible in the finite element model of the rotor bearing system, must be less than the sta- 
bility threshold. 

Lewis Research Center, 

National Aeronautics and Space Administration, 
and 

U. S. Army Air Mobility R&D Laboratory, 

Cleveland, Ohio, August 4, 1977 
505-04. 
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APPENDIX - SYMBOLS 


A time derivative of V 
a complex amplitude 
B column matrix functions of P and h 

C square matrix describing damping in finite element model 
c modal damping 

D column matrix functions of P and h 
d scalar coefficient of polynomial of hP 
E column matrix of roundoff error 
e amplitude of absolute error 

F column matrix describing forces in finite element model 
G column matrix of rate of generation of roundoff error 
g parameters used in Runge-Kutta method 
h time step 

K square matrix describing stiffnesses in finite element model 
k modal stiffness 

M diagonal matrix describing masses in finite element model 
m modal mass 

n number of mass elements in finite element model 

P partitioned square matrix used in homogeneous first-order matrix equation of 
motion 

Q partitioned column matrix describing forcing function in first-order matrix equa- 
tion of motion 

R column matrix describing displacements of shaft centerline in finite element model 
r complex elements of R 

s complex nondimensional time step 

t time 

U time derivative of Z 
u relative error growth rate 
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V time derivative of R 

w absolute error growth rate 

Z partitioned column matrix describing dependent variables in first-order matrix 
equation of motion 

e amplitude of relative error 

f damping ratio 

\ complex (damped) eigenvalue 

fjL constant used in solution of finite difference equations 

CO natural frequency 

Subscripts: 

j any mode between 1 and n 

k power of polynomial of hP 

I various terms of difference equation 
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TABLE I. - TABULATION OF dj^, FOR VARIOUS TECHNIQUES 


(a) Forward Euler 


Term of 
difference 
equation, 
1 

Power of 
polynomial 
of hP, k 

0 

1 

0 

1 

1 


(b) Backward Euler 


Term of 

Power of 

difference 

polynomial 

equation, 

of hP, k 

1 

0 

1 

-1 

0 

1 

0 

1 

0 


(c) Centered Euler 


Term of 

Power of 

difference 

polynomial 

equation. 

of hP, k 

1 

0 

1 

-1 

0 

1/2 

0 

1 

1/2 


(d) Runge-Kutta 


Term of difference 
equation, 

1 

Power of polynomial of hP, k 

0 

1 

2 

3 

4 

0 

1 

1 

1/2 

1/6 

1/24 


(e) Milne (f) Adams 


Term of 

Power of 

difference 

polynomial 

equation. 

of hP, k 

1 

0 

1 

-1 

0 

9/24 

0 

1 

19/24 

1 

0 

-5/24 

2 

0 

1/24 


Term of 
difference 
equation, 
1 

Power of 
polynomial 
of hP, k 

0 

1 

-1 

0 

1/3 

0 

0 

4/3 

1 

1 

1/3 
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Figure 1. - Model of shaft showing complex number representation of radial displacements (distance 
between shaft centerline and axis of rotation). Real and imaginary axes fixed in space. 


18 






Product of time step and natural frequency, hu 



20 



1. Report No. 

NASA TP-1092 


2. Government Accession No. 


4. Title and Subtitle 

STABILITY OF NUMERICAL INTEGRATION TECHNIQUES 
FOR TRANSIENT ROTOR DYNAMICS 

7. Author{s) 

Albert F. Kascak 

9. Performing Organization Name and Address 

NASA Lewis Research Center and 
U.S. Army Air Mobility R&D Laboratory 
Cleveland, Ohio 44135 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D. C. 20546 

15. Supplementary Notes 


3. Recipient's Catalog No. 

5. Report Date 

November 19 77 

6. Performing Organization Code 

8. Performing Organization Report No. 

E-9252 

10. Work Unit No. 

505-04 

11. Contract or Grant No. 

13 . Type of Report and Period Covered 

Technical Paper 

14. Sponsoring Agency Code 


[ 16. Abstract 

A finite element model erf a rotor bearing system was analyzed to determine the stability limits 
of the forward, backward, and centered Euler; Runge-Kutta; Milne; and Adams numerical in- 
tegration techniques. The analysis concludes that the highest frequency mode determines the 
maximum time step for a stable solution. Thus, the number of mass elements should be mini- 
mized. Increasing the damping can sometimes cause numerical instability. For a uniform 
shaft, with 10 mass elements, operating at approximately the first critical speed, the maximum 
time step for the Runge-Kutta, Milne, and Adams methods is that which corresponds to approx- 
imately 1 degree of shaft movement. This is independent of rotor dimensions. 


17. Key Words (Suggesred by AuthorlsJ 1 

Rotors 

Rotating shafts 
Numerical integration 
Numerical stability 


18. Distribution Statement 

Unclassified - unlimited 
STAR Category 37 


19. Security Classif. (of this report) 

Unclassified 


1 _ 

20. Security Classif. (of this page} 

21 . No. of Pages 

Unclassified 

21 


22. Price* 

A02 


For sale by the National Technical Information Service. Springfield. Virginia 22161 


NASA-Langley, 1977 



National Aeronautics and 
Space Administration 


THIRD-CLASS BULK RATE 


Washington, D.C. 
20546 

Official Business 

Penalty for Private Use, $300 


Postage and Fees Paid 
National Aeronautics and 
Space Administration 
NASA-451 



8 1 1U.D, 

dept of the AX8 FORCE 
AF WEAPONS LABORATORY 
ATTN: technical LIBRARY 
KIHTLAND A.FB NH 87117 


S00903DS 

(SOL) 





POSTMASTER: Undeliverable (Section 1S8 

Postal Manual) Do Not Return 



